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Abstract. We consider the simultaneous deblurring of a set of noisy images whose point spread functions are 
, different but known and spatially invariant, and the noise is Gaussian. Currently available iterative algorithms 

■ that are typically used for this type of problem are computationally expensive, which makes their application for 
very large images impractical. We present a simple extension of a classical least-squares (LS) method where the 
multi-image deblurring is efficiently reduced to a computationally efficient single-image deblurring. In particular, 
we show that it is possible to remarkably improve the ill-conditioning of the LS problem by means of stable 
operations on the corresponding normal equations, which in turn speed up the convergence rate of the iterative 
algorithms. The performance and limitations of the method are analyzed through numerical simulations. Its 

|' connection with a column weighted least-squares approach is also considered in an appendix. 

^ ■ Key words. Methods: data analysis - Methods: statistical - Techniques: Image processing 
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^ ' 1. Introduction the images taken at high galactic latitudes are expected 

•»"H . .... ..... to be almost entirely dominated by the CMB contribu- 

k> , An important problem m image processing is the simul- .... „ , . „ , 

K*% ,,, • r r i i . i r- n ^ tion in a large range of observing frequencies (e.g., see 



taneous deblurring of a set of observed images (of a fixed nTTi — , , onM , . ,, nc , n 

. . . . , . . . . . . btolvarov et al. 2002). Since the PSFs corresponding to 

obiect , where each image is degraded by a difterent spa- , , , , . , , . , . r . , 

. .. / . . . . . , nnT1 , ^ . the images obtained at the various frequencies can be quite 

tially invariant point spread function (rbb . bor example, ..... , ,, , „ , „ 

t n , . _ . ^ dissimilar, simultaneous deblurring can be useful tor mi- 

data obtained by the Large Binocular Telescope (LBT). . ,. t ™, n , , ,. r , , 

. J . r ox • proving the extraction of CMB information from the data. 
This instrument consists of two 8.4m mirrors on a com- 

mon mount with a spacing of 14. 4m between the centers . , . . .. , , . . , , . 

n 7^ — rni, „ . . , . , . Although a vast literature on the classical problem of 

(Angel ct al. 1998). i'or a given orientation of the tele- . . . . .... , , . . 

., ..j^ ~y ,. ., i , ,. , ,i smgle-imagc deblurring is available, much less has been 

scope, the diffraction- limited resolution along the center- .. . .... ... . . ,, 

.... • i , ,i r nr. n • published on multiple-image deblurring. An excellent gen- 
to-center baseline is equivalent to that of a 22.8m mir- , . . r- — f t=-*? n i " " i ° 

, ., . , ,. . ,. oral review is Bertero & Boccacci 2000b). In addition, 



ror, while the resolution along the perpendicular direc- n h ., i _ ' 

j. .-, o • a -ii , legmark 1999) provides some material more directly re- 

tion is that of a single 8m mirror. A possible way to ob- . ; ' ' ' ,. _ . , . .. 

lated to CMB studies. But one of the more serious lim- 
itations of the methods currently available is that, al- 
though able to provide satisfactory results, they are quite 
expensive in terms of computational requirements. Thus, 



tain an image with an improved and uniform spatial res- 
olution is to simultaneously deconvolve the images taken 
with different orientations of the telescope. Another ex- 
ample is the multi- frequency observations of the Cosmic , , . , • ,•„•, , , 
,,. „ , . frnmy, , , . n r , these methods may be quite dithcult to use lor the very 
Microwave Background (CMB) obtained from satellites , . , „„«.,% , r T 

i nr a ht^ts t r . i/i i m i large images (ss 10° — 10° pixels) expected from LBI 

such as PLANCK. In fact, although some other physi- ° nr . °„ T : , , . . , „ . 

and PLANCK observations. It is clear that more efficient 



cal components are present at the microwave frequencies, 



methods must be developed; the aim of this paper is to 
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In Sect.|21we formalize the problem and propose an ef- 
ficient solution in Sect. [3 Some of its possible limitations 
are considered in Sect. 0] In Sect. we study its perfor- 
mance through numerical experiments, and also compare 
it to that of other methods currently available in the lit- 
erature. Finally, Sect. [7| closes with some final comments 
and conclusions. 



It is not difficult to see that / is the solution of the normal 
equations 



J2AjA jf = J2 A J9i> 



(6) 



where Aj denotes the transpose of the block-circulant ma- 
trix A 3 . The DFT of both sides of Eq. ||BJ| provides 



2. Formalization of the problem 

Suppose one has p observed images, 9j(x,y), j = 
1, 2, . . . ,p, of a fixed two-dimensional object fo(x, y), each 
of which is degraded by a different spatially invariant blur- 
ring operator. That is, 

9j(x,y)= K j (x-x',y-y')f (x',y')dx'dy', (1) 



where, for each j, Kj(x,y) is a space invariant PSF. 

The image restoration problem of interest is to find 
an estimate f{x,y) of fo(x,y) given the observed images 
{gj(x,y)} and the known PSFs {Kj(x, y)}. 

In practice, the experimental images are discrete NxN 
arrays of pixels (for simplicity, we assume square images) 
contaminated by noise. Therefore model (TjJ has to be re- 
cast in the form 

N-l 

gj(m,n) = Kj{m — m! ,n — n')/o(m' ,n') + Wj(m,n), 

m' ,n'— 

(2) 

where Wj{m,n) is additive Gaussian white noise. For the 
moment, we assume constant standard deviations across 
images: a W] = a w . 

If the central peak (i.e., support) of each PSF is much 
smaller than the image, and if the object does not have 
structures near the image boundaries, then the convolu- 
tion product in Eq. J5J can be well approximated by cyclic 
convolution. As is well known, such an approximation is 
quite useful since it permits rewriting Eq. J2J) as 

gj(m,n) — Kj(m,n)fo(m,n) + Wj{m,n) 1 (3) 

where the symbol " " " denotes the Discrete Fourier 
Transform (DFT). 

One approach to compute an estimate of /o is to pose 
the image restoration problem as a least-squares (LS) 
problem, which is equivalent to a maximum likelihood ap- 
proach in the case of white Gaussian noise. In order to 
outline the LS approach, it is helpful to use matrix-vector 
notation: Eq. Q can be rewritten as 



(4) 



where the arrays g, f Q and Wj are obtained by column 
stacking gj (m, n) , fo (m, n) , and Wj (m, n) , respectively. Aj 
is the block-circulant matrix defined by cyclic convolution 
with the jth PSF. The LS estimate / can then be ex- 
pressed in the form: 



/ = argmin ^ \\Ajf -g^\ 



(5) 



/(m, n) \Kj(m, n)\ 2 = K*(m, n)gj(m, n) = ?9(to, n), 
j=i i=i 

(7) 

with symbol " * " denoting complex conjugation. This 
equation shows that the multi-image LS deblurring prob- 
lem (0 is equivalent to classical one-image LS deblurring 
where the "observed" image is g(m, n) — j?(m, n) and the 
PSF is K(rn,n) = Y?j=i l^i( TO i n )| 2 - ^ s is we ^ known, 
the direct LS estimate given by the inverse DFT of 



f(m,n) 



(8) 



can be very unstable. For white noise, the mean square 
error (MSE) of © is 



e[||/-/„|| 2 ] = ^£ 



f 



(9) 



The quantities Kj(m, n) can be very close to zero at high 
frequencies because the PSFs are almost band-limited in 
practice. Consequently, the MSE can be very large, which 
indicates that the LS estimate is unstable. The main task 
of any deblurring method consists of stabilizing the solu- 
tion by somehow limiting the values of such terms. Various 
methods h ave been proposed in the lite rature (e.g., see the 
review by iBertero fc BoccaccilfeOOObT) . For example, the 
Tikhonov estimate is defined as 



arg mm 



p 

E 



•mII/I 



(10) 



where /x is a positive scalar called the regularization pa- 
rameter. The additional penalty term (weighted by /i) has 
the effect of limiting the energy of /. In fact, it is possible 
to see that 



Ej K*{m, n)gj(m, n) 



(11) 



where it is evident that the parameter /i smooths out the 
influence of the smallest values of Kj(m,n). If both the 
numerator and the second term in the denominator are 
computed and stored, then for each value of fi the compu- 
tation of requires TV 2 divisions + TV 2 sums + one two- 
dimensional DFT. This means that the Tikhonov method 
is computationally efficient, but it is not always flexible 
enough to allow the implementation of constraints such 
as positivity of the solution. 
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Iterative methods provide an alternative methodol- 
ogy that exploits the semiconvergence property of iter- 
ative inversion algorithms and can be easily implemented 
with many different types of constraints on the solution. 
Two examples (which provide constrained positive solu- 
tions) are the Projected Landweber (PL) method and 
the Iterative Space Reconst ruction Algorithm (ISRA) (see 
iBertero fc Boccacci lEoOOb]) . But, although flexible, these 
methods are expensive from the computational point of 
view, especially in situations of slow convergence (for each 
iteration, a couple of two-dimensional DFTs is required). 



3. A more efficient approach 

The main limitation of the methods mentioned in the pre- 
vious section is that they are based on the normal equa- 
tions (JZJ). The resulting deblurring problem is much more 
ill-conditioned (i.e., the PSF is much broader) than im- 
plied by each of the models J2J because the coefficients 
K{m, n) in Eq. Q are squared. The most important 
consequence is that the convergence of the iterative al- 
gorithms can be remarkably slowed with obvious conse- 
quence on the computational cost. However, by multiply- 
ing both sides of Eq. Q by a function Kj (m,n), we can 
rewrite Eq. Q in the form 

K jo {m,n)f(m,n) 

_ 9k ( m > n ) + Ej^jo [Kj (™> n )/^j n )]ffi ( TO > n ) 

i + Ei^-oO^K^IVI^oK^i 2 ] 



Least Squares PSF 



(12) 



where 



K 



j (m, n) = max[ifi(m, n), K2(m, n), 



,K p (m,n)}. 



(13) 

Here the operator max[ ] extracts the element in the ar- 
ray with the largest modulus, and jo is the value of the 
corresponding index j. In the case that two or more PSFs 
whose DFTs, for a given frequency (m,n), have the same 
modulus, it is sufficient to choose one of the Kj(m, n) ac- 
cording to a rule that preserves the symmetry of the DFT 
of a real PSF. As already noted for model J7J), model (|f 2[) 
is also equivalent to a classical LS deblurring, but now the 
"observed" image is ((m,n) and the PSF is Kj (m,n). 

By construction, Kj (to, n) is the PSF whose energy 
distribution is maximally spread in the Fourier domain 
and therefore it corresponds to a PSF whose energy dis- 
tribution is maximally concentrated in the spatial domain. 
This means that C,(m,n) can be considered a sort of de- 
blurred version of i?(m, n) (see Fig. Q). In other words, 
the conditioning of system i|I2l) is much better than that 
of system @ (this is because Kj (to, n) is not obtained by 
squaring the coefficients Kj{m,n)). The important point 
is that this result has been obtained through stable op- 
erations. From these considerations, improvements in the 
computational costs are to be expected. In addition, by 
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Fig. 1. Comparison between the PSF corresponding to 
Eq. (left panel) with that corresponding to Eq. HI2JI 
(right panel) for the numerical experiments presented in 
Sect. 
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Fig. 2. Slice comparison of the two-dimensional classi- 
cal Hanning and rectangular windows with the modified 
Harming window (N w = 20). 

formulating the problem in terms of single-image deblur- 
ring, there is a gain in flexibility as standard deblurring 
algorithms can be used. 

In Appendix^ it is shown that £(to, n) is equivalently 
obtained by transforming the LS problem (J^J into a col- 
umn weighted LS problem. 

4. Two potential problems 

Although simple and potentially very effective, the proce- 
dure presented in the previous sections has two potential 
drawbacks. 

The first concerns the statistical properties of the noise 
component of £(m, n). In particular, even if the original 
noise Wj is white (i.e., a process with flat spectrum), the 
noise component of C( m J n ) nas a power-spectrum pro- 
portional to \Kj (to, n)\ 2 / Y^j \Kj(m,n)\ 2 . This is a direct 
consequence of having used the normal equations Q to 
solve the LS problem JSJ) and may have deleterious effects 
on automatic methods for stopping the iterations or es- 
timating the regularization parameters since they assume 
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Fig. 3. Original image, blurred version and blurred + 
noise version of the first image of the set (see text) and 
corresponding PSF. The images are 256 x 256 pixels. 

that noise is white. For example, in our simulations the 
generalized cross validation (GCV) method (which may be 
used to estim ate the optim al value of /i in the Tikhonov 
approach, see IVogelll2002(l consistently produced values 
that were too small. 

The second problem concerns edge effects and can be 
serious for images containing important details near their 
borders. In this case the circulant approximation used in 
model (0} is not appropriate and edge effects are to be ex- 
pected. In particular, the mean image C,(m,n) of Eq. (|12fl 
may have no relationship with the true image. We have 
to stress that this kind of problem is unavoidable and is 
shared by any other method. Its classical solution is the 
dowing of i mage a to reduce edge discontinuities. 
Vio et ai~l l)2003bj) show that a remarkable reduction 
of the edge effects in deblurring operations, that causes lit- 
tle distortion in the final result, can be obtained if, instead 
of g, we use o = g ■ h (" • " denotes element wise multi- 
plication) and h is a modification of the classical Hanning 
window 

J0.25 x a x /?, < m,n< N w ; 
0.25 x a x /J, N - N w < m, n< N; (14) 
1 otherwise. 



Fig. 4. Mean image C( TO ? n ) an( i corresponding PSF for 
the image shown in Fig. [31 deblurred version of C,{m, n), 
and deblurred version of the first image of the set (see 
text). The deblurring was done with CGLS. The images 
shown are the ones with the smallest standard deviation 
of the true residuals; for the mean image and first image of 
the set these are, respectively, 6.57 x 10 -2 and 8.67 x 10~ 2 . 



The parameters a = [1 — cos(mir / N w )} and (3 = [1 — 
cos(mr/N w )] are chosen so that the pixels in the central 
subimage are not modified and the image has continuous 
first derivatives at the edges (see Fig- EJ ■ This window ap- 
proaches the classical two-dimensional Hanning window 
as N w — * N/2 and tends to the rectangular window as 
N w — > 0. The parameter N w thus determines the filter- 
ing characteristics, in particular the frequency pass-band, 
of the window. Its "optimal" value depends on many fac- 
tors such as the noise level, the form of the PSF and the 
specific realization of the process. F or a Gaussian PSF, 
simulation results l|Vio et al. Il2003b|) show that values of 
three or four times the dispersion of the PSF provide rea- 
sonably good results. The only unavoidable consequence 
of the procedure is that a border of thickness N w has to 
be removed from the final deblurred image. 

According to these results, 6(m, n) should permit the 
use of (fl"2"|l also in situations where the circulant approx- 
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Fig. 5. Mean image ((m, n) and corresponding PSF for 
the image shown in Fig. [31 deblurred version of £(m,n), 
and deblurred version of the first image of the set (see 
text). The deblurring was done with MRNSD. The images 
shown are the ones with the smallest standard deviation 
of the true residuals; for the mean image and the first 
image of the set these are 5.92 x 10~ 2 and 8.49 x 10 -2 , 
respectively. 



imation for the blurring operator is not appropriate. We 
test this via numerical simulations. 



5. Some numerical experiments 

We present the results of some numerical simulations to 
show the flexibility and reliability of the methodology pre- 
sented above. Figs. 1-ilSI show the deblurring results for 
an extended object and a set of point-like objects. We 
have chosen these particular examples since, being non- 
smooth functions, their restoration is a difficult problem. 
Eight images are available for each object. In each case the 
PSF is a bidimcnsional Gaussian with dispersion along 
the major axis set to twelve pixels, and to four pixels 
along the minor axis. Their inclinations take equispaced 
values in the range 0° — 160°. Gaussian white noise is 
added with a standard deviation of 2% of the maximum 
value of the blurred images. Two deblurring methods have 



Fig. 6. Original image, blurred version and blurred + 
noise version of the first image of the set (see text) and 
corresponding PSF. The images are 256 x 256 pixels. 



been used: a class ical conjugate gradient method for LS 
problems (CGLS) IjBiorck Ill996j) . and a modified residual 
norm steepest descent method (MRNSD ) that can be used 
to en force nonnegativity of the solution l|Nagv k. Strakosl 
2000). For comparison, the deblurring of the first image 
of the set (that with the PSF with major axis inclined at 
0°) is also shown. In every case the improvement due to 
the use of the eight images is evident. 

The two examples just presented do not suffer edge 
effects. The objects are in the center of the images and 
their borders show only "dark sky" . A more difficult ex- 
ample is provided by images with details near their bor- 
ders since in this case the circulant approximation for 
the blurring operator is incorrect. In order to show the 
results obtainable in this case, we consider observations 
of the Cosmic M icrowave Background (CMB) (e.g., see 
IVio et al. 1120034 - Although the PSFs considered in this 
work are different from those typical of CMB observations, 
we have chosen this example because the CMB signal is 
expected to be the realization of a stationary stochastic 
process covering the whole sky. It represents one of the 
most unfavourable situations in the deblurring of astro- 
nomical images. Figs. l9ll21 show that, despite the difficulty 



6 



R. Vio, J. Nagy, L. Tenorio, & W. Wamsteker: Multiple image deblurring 
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Fig. 7. Mean image ((m, n) and corresponding PSF for 
the image shown in Fig. deblurred version of £(m,n), 
and deblurred version of the first image of the set (see 
text). The deblurring was done with CGLS. The images 
shown are the ones with the smallest standard deviation 
of the true residuals; for the mean image and the first 
image of the set these are 6.09 x 10~ 3 and 6.17 x 10 -3 , 
respectively. 



Fig. 8. Mean image C( TO ? n ) an( i corresponding PSF for 
the image shown in Fig. |3 deblurred version of C( m j n )i 
and deblurred version of the first image of the set (see 
text). The deblurring was done with MRNSD (1000 iter- 
ations). The images shown are the ones with the smallest 
standard deviation of the true residuals; for the mean im- 
age and the first image of the set these are 5.52 x 10~ 3 
and 6.02 x 10 -3 , respectively. 



of the problem, the application of the windowing proce- 
dure presented in Sect. 01 can provide good results. The 
experiments have been carried out under the same con- 
ditions as the previous ones, with the difference that two 
different levels of noise have been tested (2% and 10% of 
the dispersion of the blurred signal) , and the method used 
in the deblurring is the classical Tikhonov app roach. The 
CMB m aps have been simulated as described in I Vio et al. I 
2003af) . Again the results are reasonably good. 



Our examples have shown that the use of C( to j n) cou- 
pled with standard deblurring algorithms can remarkably 
improve the results. However, it is still necessary to check 
the effective gain in computational cost. Figs. Upland ITT1 
show that the convergence rate of the MRNSD algorithm 
applied to C( m ; n ) is faster than that of the PL or ISRA 
methods applied to i?(m, n) (the cost per iteration of each 
method is essentially the same) . We remark that MRNSD 
did not perform very well when applied to the more ill- 



conditioned problem involving $(m, n). This is not too 
surprising since MRNSD is essentially a steepest descent 
method, which is known to have poorer convergence prop- 
er ties for more ill-condition ed problems. It has been shown 
in lNaev fe StrakoFI J2000) that preconditioning can dra- 
matically increase the convergence rate for MRNSD, but 
it is a nontrivial process to choose an appropriate precon- 
diti oner, especially for seve rely ill-conditioned problems; 
see lNagy fc Strakoil 1 2000|) for further details. 

Furthermore, to check that this result is not due to 
the kind of algorithm used, Figs. ED and 1161 compare the 
convergence rate when the PL and ISRA methods are ap- 
plied to both C(to, n ) and "&{m, n). Again, it is evident that 
the use C( m i n ) can improve the convergence rate of the 
algorithms. 
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Ordinal Image PSF of tte First I mags Wlean Imaga Maan PSF 




Fig. 9. Original image, blurred version and blurred + 
noise version of the first image of the set (see text) and 
corresponding PSF. The images are 340 x 340 pixels. Noise 
ofS/N = 2. 

6. Extension to images with different noise levels 

So far we have assumed that the noise level is the same 
in all the images (cr Wj = cr w ), a condition that allows the 
combination of different images to provide the most inter- 
esting results. Moreover, this requirement is often satisfied 
in practice. For example, with LBT the images taken at 
different orientations of the telescope are expected to be 
characterized by the same noise level. 

The extension of the method to images with differ- 
ent noise levels is straightforward. If the different a Wi are 
known, one may account for the different variabilities by 
changing equation (0 to 

/ = argmin £ || (A,f - 9j )/a Wj f, (15) 

3=1 

or, alternatively, 

p 

f = argmin £ \\Ajf - Qj f, (16) 

3=1 

where A.j and represent, respectively, the matrix Aj 
and array with their elements divided by a Wj . The rest 



Fig. 10. Mean image (,{m : n) and corresponding PSF 
for the image shown in Fig. ^\ deblurred version of 
C(m,n), and deblurred version of the first image of the 
set (see text). The deblurring was done with the classical 
Tikhonov method and the windowing procedure explained 
in the text (N w = 30 pixels). The images shown are the 
ones with the smallest standard deviation of the true resid- 
uals that for the mean image and the first image of the 
set are 0.40 and 0.52, respectively. 

of the procedure remains as described in Sect. For ex- 
ample, Eq. ifHfll becomes 

lCj (to, n) = max [ /Ci (m, n),fC2(m, n), . . . , K, p (m, n) ], 

(17) 

where ICj (to, n) = Kj (m, n) / a Wj . 
7. Concluding remarks 

We have considered the problem of simultaneously deblur- 
ring a set of images of a fixed object blurred by different 
space-invariant PSFs. Although currently available meth- 
ods seem to provide good results, they are comptutation- 
ally expensive. We have developed a method, based on a 
LS approach, to efficiently transform a multi-image de- 
blurring problem into a single-image one. This approach 
provides substantial savings in computational require- 
ments and can be implemented using standard currently 
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Ordinal Image PSF of tte First I mags Wlean Imaga Maan PSF 




Fig. 11. Original image, blurred version and blurred + 
noise version of the first image of the set (see text) and 
corresponding PSF. The images are 340 x 340 pixels. Noise 
of S/N = 10. 

available numerical algorithms. These conclusions are con- 
firmed by our numerical experiments. 

But despite these encouraging results, some questions 
are still open. In particular, regularization parameter se- 
lection methods (such as GCV for the Tikhonov approach) 
have to be extended to deal with the correlated noise com- 
ponent of the mean image. This is especially important for 
the deblurring of random fields where the statistical prop- 
erties of the field are more important than a particular 
realization. 

Another question is the generalization of the proposed 
method to non-Gaussian noise. For example, for Poisson 
noise good results have been obtained with the maximum- 
likelihood estimate 

p N-l 

f = arg max ^ 

j—1 m,n—0 

[g j (m,n)ln(A j f)(m,n) - (Aj /)(m,n)], (18) 

l|Correia fc Richichi lliooot iBertero fc Boccacci lE000a|) . 

Acknowledgements. We thank Prof. M. Bertero (Universita di 
Genova) for useful discussions. 



Fig. 12. Mean image (,{m : n) and corresponding PSF 
for the image shown in Fig. 1111 deblurred version of 
£(m,n), and deblurred version of the first image of the 
set (see text). The deblurring was done with the classical 
Tikhonov method and the windowing procedure explained 
in the text (N w = 30 pixels). The images shown are the 
ones with the smallest standard deviation of the true resid- 
uals that for the mean image and the first image of the 
set are 0.35 and 0.48, respectively. 

Appendix A: Multi-Frame Deblurring and 
Weighted Least-Squares 

In this appendix the multi-frame deblurring problem pre- 
sented in the previous sections is shown to be related to 
(column) weighted LS. 

In the multi-frame deblurring problem, we consider the 
LS problem (we omit regularization, but this can be easily 
incorporated into the discussion): 
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Assuming each Aj is block circulant, and can be diago- 
nalized by the unitary Fourier transform matrix, then this 
problem can be transformed to the equivalent LS problem: 
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Fig. 13. Deblurred images, corresponding to the image 
shown in Fig. [31 obtained through the MRNSD method 
applied to ((m, n) , and the PL, and ISRA methods ap- 
plied to $(m, n). The images shown are the ones with the 
smallest standard deviation of the true residuals that are 
5.92 x KT 2 , 5.39 x 1(T 2 , and 5.63 x 1(T 2 respectively. 
The convergence curves for each method is also shown: 
MRSND converges after 346 iterations, whereas both PL 
and ISRA do not converge after 1500 iterations. 
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(A.2) 



where Aj = diag(A^, . . . , A^ ; ) is a diagonal matrix 
containing the eigenvalues of the matrix Aj (i.e., this is 
just the DFT of the PSF, or in the notation of the main 
text Kj). 

Now, following the arguments in Sect. let Si — 

max {Ap^} , i = 1, 2, . . . , N. Define the diagonal matrix: 

i<j<p 



vO') 



Fig. 14. Deblurred images, corresponding to the image 
shown in Fig. [fj] obtained through the MRNSD method 
applied to £(m, n) , and the PL, and ISRA methods ap- 
plied to t?(m, n). The images shown are the ones with the 
smallest standard deviation of the true residuals that are 
5.52 x 10~ 3 , 6.05 x 10~ 3 , and 5.69 x 10~ 2 respectively. The 
convergence rate for each method is also shown: none of 
the methods converged before 1000 iterations. 

and consider the (column) weighted LS problem: 



(A.4) 
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We will write this LS problem as: 

min || DA/- 9 II2, 

where 



D 



AiA 

A 2 A _1 



ApA -1 



(A.5) 



(A.6) 



diag((Si,<5 2 , ...,S N ). 



(A.3) 



The normal equations version of this LS problem is: 

A* D* DAf = A*D*g , (A.7) 
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Fig. 15. Convergence rate, corresponding to the image 
shown in Fig. [31 for the PS and ISRA methods applied 
to both C(m,n) and i d(m,n). 
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Fig. 16. Convergence rate, corresponding to the image 
shown in Fig. for the PS and ISRA methods applied 
to both £(m,n) and i d(m,n). 



or, equivalently, 



D*DAf = D*g, 



(A. 



where D* is the complex conjugate transpose of D. 

Now, observe that D* D is perfectly well conditioned, 
and thus can be easily inverted. To see this, note that for 
all i, there exists jo = jo(i) such that 



and \X^\>\\^\, j^jo. 



Si = Ap o) 
The ith diagonal entry of D*D is 
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(A.9) 



(A.10) 



and thus each diagonal entry of D*D can be bounded by 



Ma« I 2 
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(A.ll) 



Thus D*D is perfectly well conditioned. 

Since D*D is perfectly well conditioned, there is no 
danger in transforming the normal equations given in 
equation HA.8|) to 



Af = (D*D) 1 D*g. 



(A.12) 



which is equivalent to Eq. (|12fl . 

A vast literature is a vailable on column w eighte d 
LS; see, for example, iGolub k. Van LoarTI ljl996j) . 
lLawson k Hanson I l)l995|) . and lBiorck I l|l996|) . 
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